	subroutine dsetup( n, V, pd, nr, normV, full, band, b, job )
	integer           n, job, band, pd(*), nr(*)
	double precision  V(*), b(*)
	double precision  normV
	logical           full

c       ... Local Variables ...
	integer          iv
	double precision nr_sum, nr_mean, nr_std, nr_var
	double precision normcol

*------------------------------------------------------------------------------
*	DSETUP sets up the array V corresponding to the matrix A. The matrix
*              is stored in a packed fishbone fashion as described below.
*              It also sets up the rigt hand side.
*
*------------------------------------------------------------------------------
*
*	Input:
*
*	n :   INTEGER
*             size of the problem
*
*       full: LOGICAL
*             if .TRUE.  =>  generate a full n-by-n matrix.
*             if .FALSE. =>  generate a sparse matrix.
*
*       band: INTEGER
*             Indicates what sort of skyline to generate and the size
*             of the band compared to the dimension of the problem.
*
*             band = 10  =>  generate a uniform banded matrix, where the
*                            nonzero elements are 10% of the dimension of
*                            the problem.
*
*             band <= 0   
*                 or     =>  generate random band.
*             band > 100
*
*             If full = .TRUE., band is ignored.
*
*	job : INTEGER
*             Indicates the data to generate.
*             job = 1  =>  generate matrix V only.
*             job = 2  =>  generate rhs b only.
*             job = 3  =>  generate both matrix V and rhs b.
*
*
*       Output:
*
*       b : the right hand side
*
*	V : array where the elements are stored in a packed fashion following
*	    a fishbone design. That way lines below the diagonal, and columns 
*	    above the diagonal are stored with stride 1.
*	    The matrix ma has a symmetrical kind of banded-like shape. But
*	    the values of the matrix are not symmetrical.
*
*       pd : an integer array containing pointers to the diagonal elements
*            of the matrix A in the array V.
*	   
*       nr : an integer array containing the number of nonzero elements in
*            a row or column up to the diagonal.
*	   
*
* V                                A
* ------------------------------   ------------------------------------
* |  1 | 3 | 7 | * | * | * | * |   | 11 | 12 | 13 |  * |  * |  * |  * |
* |----|   |   |   |   |   |   |   |----|    |    |    |    |    |    |
* |  2 | 4 | 8 |12 |18 | * | * |   | 21 | 22 | 23 | 24 | 25 |  * |  * |
* |--------|   |   |   |   |   |   |---------|    |    |    |    |    |
* |  5   6 | 9 |13 |19 | * |31 |   | 31 | 32 | 33 | 34 | 35 |  * | 37 |
* |------------|   |   |   |   |   |--------------|    |    |    |    |
* |  *  10  11 |14 |20 |24 |32 |   |  *   42   43 | 44 | 45 | 46 | 47 |
* |----------------|   |   |   |   |-------------------|    |    |    |
* |  *  15  16  17 |21 |25 |33 |   |  *   52   53   54 | 55 | 56 | 57 |
* |--------------------|   |   |   |------------------------|    |    |
* |  *   *   *  22  23 |26 |34 |   |  *    *    *   64   65 | 66 | 67 |
* |------------------------|   |   |-----------------------------|    |
* |  *   *  27  28  29  30 |35 |   |  *    *   73   74   75   76 | 77 |
* ------------------------------   ------------------------------------
*
*
*	pd : array of diagonal indexes
*		(e.g. 1, 4, 9, 14, 21, 26, 35)
*
*
*       For a full matrix:
*
*		------------------------------------
* 		|  1 |  3 |  7 | 13 | 21 | 31 | 43 | 
* 		| 11 | 12 | 13 | 14 | 15 | 16 | 17 |
* 		|---------|    |    |    |    |    |
* 		|  2 |  4 |  8 | 14 | 22 | 32 | 44 |
* 		| 21 | 22 | 23 | 24 | 25 | 26 | 27 |
* 		|---------------    |    |    |    |
* 		|  5    6 |  9 | 15 | 23 | 33 | 45 |
* 		| 31   32 | 33 | 34 | 35 | 36 | 37 |
* 		|--------------------    |    |    |
* 		| 10   11   12 | 16 | 24 | 34 | 46 |
* 		| 41   42   43 | 44 | 45 | 46 | 47 |
* 		|-------------------------    |    |
* 		| 17   18   19   20 | 25 | 35 | 47 |
* 		| 51   52   53   54 | 55 | 56 | 57 |
* 		|------------------------------    |
* 		| 26   27   28   29   30 | 36 | 48 |
* 		| 61   62   63   64   65 | 66 | 67 |
* 		|----------------------------------|
* 		| 37   38   39   40   41   42 | 49 |
* 		| 71   72   73   74   75   76 | 77 |
* 		------------------------------------
*
*	pd = ( 1, 4, 9, 16, 25, 36, 49 )
*
*------------------------------------------------------------------------------


c   ... Generate PD ...

	init = 1325

c   ... Full Matrix ...
	if ( full ) then
	   do i = 1, n
	      pd(i) = i*i
	   end do

        else
	   pd(1) = 1

           if ( band .LE. 0  .OR.  band .GT. 100 ) then

c         ... Random bandwidth ...
	      do i = 2, n
	         pd(i) = pd(i-1) + 2*mod( 100*ran(init), i ) + 1
              end do

           else

c         ... Fixed bandwidth ...
              nband = n*band / 100
              do i = 2, n
                 pd(i) = pd(i-1) + min( nband + 1, 2*(i-1) + 1 )
              end do
           end if

	end if


c   ... Calculate NR (number of nonzero elements in a row up to      ...
c   ... the diagonal), its average value and the standard deviation. ...

	nr_sum = 0
	nr(1)  = 0
	do i = 2, n
	   nr(i)  = ( pd(i) - pd(i-1) ) / 2
	   nr_sum = nr_sum + nr(i)
	end do
	nr_mean = nr_sum / float(n-1)

	nr_sum = 0
	do i = 1, n
	   nr_sum = nr_sum + ( nr(i) - nr_mean )**2
        end do
	nr_var = nr_sum / float(n-1)
	nr_std = sqrt( nr_var )


c   ... Generate  a matrix that will not require pivoting ...

	init = mod( 3125*init, 65536 )
	V( 1 ) = ( init - 32768 ) / 16384
	do j = 2, n
	   do i = pd(j-1)+1, pd(j)
	      init = mod( 3125*init, 65536 )
	      V( i ) = ( init - 32768 ) / 16384
	   end do
	end do
	do j = 1, n
	   normcol = 0.d0
	   do i = pd(j) - nr(j), pd(j)
	      normcol = normcol + abs( V( i ) )
	   end do
	   do i = j+1, n
	      if( nr(i) .GE. i-j ) 
     &            normcol = normcol + abs( V( pd(i) - nr(i) -(i-j) ) )
	   end do
	   V( pd( j ) ) = V( pd( j ) ) + normcol + 1.d0
	end do


c  ...  Calculate the norm of A. ...

	normV = max( V( 1 ), 0.d0 )
	do j = 2, n
	   do i = pd( j-1 ) + 1, pd( j )
	      normV  = max( V( i ), normV )
	   end do
	end do


c   ... Generate the right hand side ...

	do i = 1, n
	   b(i) = 0.d0
	end do

	do j = 1, n
	   do i = pd(j) - nr(j), pd(j)
	      b( j - pd(j) + i ) = b( j - pd(j) + i ) + V(i)
	   end do
	   do i = j+1, n
	      if( nr(i) .GE. i-j )
     &	          b(i) = b(i) + V( pd(i) - nr(i) -(i-j) )
	   end do
	end do


	return
	end
